function srf = vmr_DBReco(hfile, bcol, scol)
% VMR::DBReco  - direct border reconstruction
%
% FORMAT:       srf = vmr.DBReco([bcol, scol])
%
% Input fields:
%
%       bcol        border color (default 235)
%       scol        segmentation color (default 240)
%
% Output fields:
%
%       srf         reconstructed surface

% Version:  v0.7c
% Build:    7101317
% Date:     Oct-13 2007, 5:18 PM CEST
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
if nargin < 2 || ...
    numel(bcol) ~= 1 || ...
    ~isnumeric(bcol) || ...
    isnan(bcol) || ...
    bcol < 226 || ...
    bcol > 255 || ...
    bcol ~= fix(bcol)
    bcol = uint8(235);
else
    bcol = uint8(bcol);
end
if nargin < 3 || ...
    numel(scol) ~= 1 || ...
    ~isnumeric(scol) || ...
    isnan(scol) || ...
    scol < 226 || ...
    scol > 255 || ...
    scol ~= fix(scol)
    scol = uint8(240);
else
    scol = uint8(scol);
end

% get VMR content
bc = bvqxfile_getcont(hfile.L);
vd = bc.VMRData(:, :, :);
if ~isa(bc.VMRData, 'uint8')
    error( ...
        'BVQXfile:BadCall', ...
        'DBReco method only valid for 8-bit VMR data.' ...
    );
end
vs = size(vd);
vsp = vs + 1;

% re-prepare
vd(vd == bcol) = scol;

% create boolean matrix with non matching color
vdn = (vd ~= scol);

% find border by notting array with dilation operator
vd(~vdn & conv3d(vdn, 2)) = bcol;
[vdn, vds] = flexmask(vd, double(bcol), 0, 1, vs + 2);

% make sure to get biggest chunk
[vdi{1:3}] = ind2sub(size(vdn), find(vdn(:)));
vdf = floodfill3(vdn, vdi{1}(1), vdi{2}(1), vdi{3}(1), 'edge');
while (sum(vdf(:)) < (0.5 * vds))
    vdn(vdf) = false;
    if sum(vdn(:)) == 0
        error( ...
            'BVQXfile:InternalError', ...
            'No more marked voxels.' ...
        );
    end
    [vdi{1:3}] = ind2sub(size(vdn), find(vdn(:)));
    vdf = floodfill3(vdn, vdi{1}(1), vdi{2}(1), vdi{3}(1), 'edge');
end
vdn(~vdf) = false;

% create coordinate list array
clst = uint32([]);
clst(vsp(1), vsp(2), vsp(3)) = 0;

% find list of voxels that are marked
[cx, cy, cz] = ind2sub(size(vdn), find(vdn));
vdn = flexmask(vd == bcol | vd == scol, 0, 2, 0, vs + 2);
cxyz = [cx(:), cy(:), cz(:)];

% at most, each voxel can produce 12 triangles and 8 points
tlist = uint32([]);
tlist(12 * numel(cx), 3) = 0;

% initialize counter
cc = 1;
tc = 1;
pl = uint32(zeros(1, 4));

% iterate over found voxels
for vc = 1:numel(cx)
    
    % for fast access get cube with data
    crdf = cxyz(vc, :) - 1;
    crdt = crdf + 2;
    cu = vdn(crdf(1):crdt(1), crdf(2):crdt(2), crdf(3):crdt(3));
    
    % check faces, ~ means no border!
    % face in direction of A (-x)
    if ~cu(1, 2, 2)
        
        % get center of two triangle faces
        crd = cxyz(vc, :) - [1, 0.5, 0.5];
        
        % get four coordinates
        p = [ ...
            crd + [0, -0.5, -0.5]; ...
            crd + [0,  0.5, -0.5]; ...
            crd + [0, -0.5,  0.5]; ...
            crd + [0,  0.5,  0.5]];
        
        % check whether a point exists or add
        for pc = 1:4
            ppnum = clst(p(pc, 1), p(pc, 2), p(pc, 3));
            if ppnum > 0
                pl(pc) = ppnum;
            else
                pl(pc) = cc;
                clst(p(pc, 1), p(pc, 2), p(pc, 3)) = cc;
                cc = cc + 1;
            end
        end
        
        % check direction where this point goes to
        
        % create two triangles
        tlist(tc, 1:3) = pl([1, 2, 4]);
        tlist(tc + 1, 1:3) = pl([4, 3, 1]);
        tc = tc + 2;
    end
    
    % face in direction of P (+x)
    if ~cu(3, 2, 2)
        
        % get center of two triangle faces
        crd = cxyz(vc, :) - [0, 0.5, 0.5];
        
        % get four coordinates
        p = [ ...
            crd + [0, -0.5, -0.5]; ...
            crd + [0,  0.5, -0.5]; ...
            crd + [0, -0.5,  0.5]; ...
            crd + [0,  0.5,  0.5]];
        
        % check whether a point exists or add
        for pc = 1:4
            ppnum = clst(p(pc, 1), p(pc, 2), p(pc, 3));
            if ppnum > 0
                pl(pc) = ppnum;
            else
                pl(pc) = cc;
                clst(p(pc, 1), p(pc, 2), p(pc, 3)) = cc;
                cc = cc + 1;
            end
        end
        
        % create two triangles
        tlist(tc, 1:3) = pl([1, 4, 2]);
        tlist(tc + 1, 1:3) = pl([4, 1, 3]);
        tc = tc + 2;
    end

    % face in direction of S (-y)
    if ~cu(2, 1, 2)
        
        % get center of two triangle faces
        crd = cxyz(vc, :) - [0.5, 1, 0.5];
        
        % get four coordinates
        p = [ ...
            crd + [-0.5, 0, -0.5]; ...
            crd + [-0.5, 0,  0.5]; ...
            crd + [ 0.5, 0, -0.5]; ...
            crd + [ 0.5, 0,  0.5]];
        
        % check whether a point exists or add
        for pc = 1:4
            ppnum = clst(p(pc, 1), p(pc, 2), p(pc, 3));
            if ppnum > 0
                pl(pc) = ppnum;
            else
                pl(pc) = cc;
                clst(p(pc, 1), p(pc, 2), p(pc, 3)) = cc;
                cc = cc + 1;
            end
        end
        
        % create two triangles
        tlist(tc, 1:3) = pl([1, 4, 3]);
        tlist(tc + 1, 1:3) = pl([4, 1, 2]);
        tc = tc + 2;
    end
    
    % face in direction of I (+y)
    if ~cu(2, 3, 2)
        
        % get center of two triangle faces
        crd = cxyz(vc, :) - [0.5, 0, 0.5];
        
        % get four coordinates
        p = [ ...
            crd + [-0.5, 0, -0.5]; ...
            crd + [-0.5, 0,  0.5]; ...
            crd + [ 0.5, 0, -0.5]; ...
            crd + [ 0.5, 0,  0.5]];
        
        % check whether a point exists or add
        for pc = 1:4
            ppnum = clst(p(pc, 1), p(pc, 2), p(pc, 3));
            if ppnum > 0
                pl(pc) = ppnum;
            else
                pl(pc) = cc;
                clst(p(pc, 1), p(pc, 2), p(pc, 3)) = cc;
                cc = cc + 1;
            end
        end
        
        % create two triangles
        tlist(tc, 1:3) = pl([1, 3, 2]);
        tlist(tc + 1, 1:3) = pl([3, 4, 2]);
        tc = tc + 2;
    end
    
    % face in direction of R (-z)
    if ~cu(2, 2, 1)
        
        % get center of two triangle faces
        crd = cxyz(vc, :) - [0.5, 0.5, 1];
        
        % get four coordinates
        p = [ ...
            crd + [-0.5, -0.5, 0]; ...
            crd + [-0.5,  0.5, 0]; ...
            crd + [ 0.5, -0.5, 0]; ...
            crd + [ 0.5,  0.5, 0]];
        
        % check whether a point exists or add
        for pc = 1:4
            ppnum = clst(p(pc, 1), p(pc, 2), p(pc, 3));
            if ppnum > 0
                pl(pc) = ppnum;
            else
                pl(pc) = cc;
                clst(p(pc, 1), p(pc, 2), p(pc, 3)) = cc;
                cc = cc + 1;
            end
        end
        
        % create two triangles
        tlist(tc, 1:3) = pl([1, 4, 2]);
        tlist(tc + 1, 1:3) = pl([4, 1, 3]);
        tc = tc + 2;
    end
    
    % face in direction of L (+z)
    if ~cu(2, 2, 3)
        
        % get center of two triangle faces
        crd = cxyz(vc, :) - [0.5, 0.5, 0];
        
        % get four coordinates
        p = [ ...
            crd + [-0.5, -0.5, 0]; ...
            crd + [-0.5,  0.5, 0]; ...
            crd + [ 0.5, -0.5, 0]; ...
            crd + [ 0.5,  0.5, 0]];
        
        % check whether a point exists or add
        for pc = 1:4
            ppnum = clst(p(pc, 1), p(pc, 2), p(pc, 3));
            if ppnum > 0
                pl(pc) = ppnum;
            else
                pl(pc) = cc;
                clst(p(pc, 1), p(pc, 2), p(pc, 3)) = cc;
                cc = cc + 1;
            end
        end
        
        % create two triangles
        tlist(tc, 1:3) = pl([3, 2, 4]);
        tlist(tc + 1, 1:3) = pl([2, 3, 1]);
        tc = tc + 2;
    end
end

% set pointers correctly
tc = tc - 1;
cc = cc - 1;

% cut lists
tlist = double(tlist(1:tc, :));
cfnd = find(clst);
[cidx{1:2}] = sort(clst(cfnd));
cfnd = cfnd(cidx{2});
[cx, cy, cz] = ind2sub(vsp, cfnd);
if bc.VoxResX >= 0.5
    cxyz = [ ...
        bc.VoxResX * (cx(:) + bc.OffsetX - 1.5), ...
        bc.VoxResY * (cy(:) + bc.OffsetY - 1.5), ...
        bc.VoxResZ * (cz(:) + bc.OffsetZ - 1.5)];
else
    cxyz = [ ...
        cx(:) + bc.OffsetX - 1.5, ...
        cy(:) + bc.OffsetY - 1.5, ...
        cz(:) + bc.OffsetZ - 1.5];
end

% get neighbors
[nei, bn] = mesh_trianglestoneighbors(cc, tlist);

% warn if bad neighborhood
fixc = 0;
if ~isempty(bn)
    warning( ...
        'BVQXfile:InternalError', ...
        'Error with neighborhood relation, attempting to fix.' ...
    );
end
while ~isempty(bn)
    fixc = fixc + 1;
    if fixc > 3
        error( ...
            'BVQXfile:InternalError', ...
            'Neighborhood problem unfixable.' ...
        );
    end
    cxyz = [cxyz; cxyz(bn, :)];
    for bnc = 1:numel(bn)
        cc = cc + 1;
        bni = bn(bnc);
        cnei = nei{bni, 2};
        atr = find(any(tlist == bni, 2));
        at = tlist(atr, :);
        for atc = 1:size(at, 1)
            if isempty(intersect(cnei, at(atc, :)))
                tlist(atr(atc), at(atc, :) == bni) = cc;
            end
        end
    end
    [nei, bn] = mesh_trianglestoneighbors(cc, tlist);
end

% make sure that all vertices are connected (only one surface)
neidone = false(1, cc);
neicons = false(1, cc);
while sum(neidone) < (0.5 * cc)
    neir = nei(:, 2);
    neicons(1) = true;
    while any(neicons)
        neidone(neicons) = true;
        consnei = neir(neicons);
        neicons(:) = false;
        for c = 1:numel(consnei)
           neicons(consnei{c}) = true;
        end
        neicons(neidone) = false;
    end
    if sum(neidone) >= (0.5 * cc)
        rmi = find(~neidone);
    else
        rmi = find(neidone);
    end
    if ~isempty(rmi)
        tii = zeros(cc, 1);
        cc = cc - numel(rmi);
        cxyz(rmi, :) = [];
        for coli = 1:3
            isi2 = 1;
            while ~isempty(isi2)
                [isr, isi1, isi2] = intersect(rmi, tlist(:, coli));
                tlist(isi2, :) = [];
            end
        end
        tiu = unique(tlist(:));
        tii(unique(tlist(:))) = 1:numel(tiu);
        tlist = tii(tlist);
        [nei, bn] = mesh_trianglestoneighbors(cc, tlist);
        if ~isempty(bn)
            error( ...
                'BVQXfile:InternalError', ...
                'Removing separate surface corrupted triangle list.' ...
            );
        end
    end
end
tc = size(tlist, 1);

% compute Euler characteristic
eul = cc - 0.5 * tc;
if eul ~= 2
    warning( ...
        'BVQXfile:BadObject', ...
        'Probably %d handles/bridges in surface remaining.', ...
        round((2 - eul) / 2) ...
    );
end

% remove neighbors that are not diagonals
for vc = 1:cc
    n = nei{vc, 2};
    ncrd = (abs(sum(cxyz(vc * ones(1, numel(n)), :) - cxyz(n, :), 2)) == 1);
    nei{vc, 1} = sum(ncrd);
    nei{vc, 2}(~ncrd) = [];
end

% build surface
srf = emptysrf;
srfc = bvqxfile_getcont(srf.L);
srfc.NrOfVertices = cc;
srfc.NrOfTriangles = tc;
srfc.VertexCoordinate = cxyz;
srfc.VertexColor = zeros(cc, 4);
srfc.Neighbors = nei;
srfc.TriangleVertex = tlist;

% put into storage array
bvqxfile_setcont(srf.L, srfc);

% calculate normals !
srf_RecalcNormals(srf);
